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The effect of a step wise change in the pihar density on the dynamics of droplets is investigated via 
three-dimensional lattice Boltzmann simulations. For the same pillar density gradient but different 
pillar arrangements, both motion over the gradient zone as well as complete arrest are observed. 
In the moving case, the droplet velocity scales approximately linearly with the texture gradient. 
A simple model is provided reproducing the observed linear behavior. The model also predicts a 
linear dependence of droplet velocity on surface tension. This prediction is clearly confirmed via 
our computer simulations for a wide range of surface tensions. 



I. INTRODUCTION 

Due to its occurrence in a wide range of natural phe- 
nomena and its fundamental importance for surface en- 
gineering [l|, m , the behavior of liquid drops on solid sur- 
faces is an active field of research [3]. Individual droplets 
can serve as ideal chemical reactors [Jl , carriers of infor- 
mation ^ or in ink-jet printers as well as surface prepa- 
ration prior to painting or coating. 

Although there exist a few solids that are molecu- 
larly fiat (e.g. mica), most of solids are rough on the 
micro scale [6] so that the^ well known Young's law, 
cos^Y = (c^sv — crsL)/crLV [7], must be modified in or- 
der to take account of surface roughness. Here, is the 
contact angle on a fiat substrate and (Jlv, ctsl and <Jsv 
are the liquid- vapor, solid-liquid and solid- vapor specific 
surface free energies, respectively. The simplest and most 
popular modification of the Young's law for rough sur- 
faces dates back to the works of Wenzel [8] and Cassie 
and Baxter 0, where the effect of roughness on wetting 
is assumed to be a mere change of the average surface 
areas involved in the problem. 

Assuming that the liquid completely penetrates into 
the roughness grooves (collapsed state), Wenzel obtained 
cos 0^ = r cos 9y for the apparent contact angle 0^ (the 
roughness factor r is the real solid area within a square 
of unit length). Cassie and Baxter, on the other hand, 
considered the case of a droplet pending on the top of 
roughness tips (suspended state) and obtained 

cos i9£ =^cos6>Y-(l-^), (1) 

where the roughness density (j) gives the fraction of the 
droplet 's base area, which is in contact with the solid. 
It is important to realize that both the Wenzel and the 
Cassie-Baxter equations do not explicitly take account 
of three phase contact line structure. This shortcoming 
may, however, be neglected as long as the contact area 
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refiects the structure and energetics of the three phase 
contact line jlQi] . 

The suspended state is often separated from the Wen- 
zel state by a finite free energy barrier, which depends 
both on the droplet size and roughness characteristics 
pj] , [12I , \\^ . On a rough hydrophobic substrate, the ap- 
parent contact angle of a droplet in the suspended state 
is typically higher than in the collapsed state [H, [13, [lH . 
Furthermore, the contact angle hysteresis significantly in- 
creases when a suspended droplet undergoes a transition 
to the Wenzel state 0, Q • Indicative of stronger pinning 
of the three phase contact line, this feature refiects 
itself in a sticky behavior of liquid drops in the collapsed 
state [2I, as compared to their high mobility in the 
suspended state. 

While droplet behavior on homogeneous roughness has 
widely been investigated in the literature, only few works 
exist dealing with the case of inhomogeneous topography 
[H, m, [2O, dll, [22[ . Indeed, experimental observation of a 
roughness gradient induced spontaneous motion is not an 
easy task [21|. The authors of [20, 21], for example, resort 
to shaking vertically the substrate in order to overcome 
pinning forces. 

Recently, a two-dimensional theoretical model is pro- 
posed aiming at a study of the present topic [22] . To the 
best of our knowledge, however, computer simulations 
of the problem are lacking so far. The present work is 
aimed at filling this gap. We provide first direct numer- 
ical evidence for spontaneous droplet motion actuated 
by a gradient of pillar density. Furthermore, we investi- 
gate the influence of specific distribution/arrangement of 
roughness elements (pillars in our case) on the behavior 
of the droplet. An important observation is that, de- 
pending on the specific arrangement of the pillars, both 
complete arrest and motion over the entire gradient zone 
can be observed for the same gradient of pillar density. 
This underlines the importance of the topography design 
for achieving high mobility drops. 

For the case of mobile drops, we provide a simple model 
for the dependence of drop velocity on pillar density dif- 
ference A<j) = (/)Right — ^Left- The model accounts for 
the observed linear dependence and predicts further that 
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the velocity should scale linearly with the liquid-vapor 
surface tension. This prediction is also in line with our 
computer simulations for the studied range of parame- 
ters. 



II. NUMERICAL MODEL 

We employ a free energy based two phase lattice Boltz- 
mann (LB) method, first proposed by Swift [23^. After 
establishing the Galilean invariance [24] , it was developed 
further (25|, l27| in order to take the wetting effect of the 
solid substrate into account. Since then, the approach 
has been used to study e.g. stability and dynamics of 
droplets on topographically patterned hydrophobic sub- 
strates [U [281] , effect of chemical surface patterning 
on droplet dynamics [H, [2^ as well as chemical gradient 
induced separation of emulsions [sO] • A detailed descrip- 
tion of the method can be found in [25l,[2?]. For the sake 
of completeness, however, we present a short overview of 
the method. Equilibrium properties of the present model 
can be obtained from a free energy functional 

*=/" [Mpiv)) + ^{dap{v)f) dr^ + [ ^p^ds. (2) 
J V ^ s 

In Eq. ([2]), ?/^b is the bulk free energy density of the fluid, 
V denotes the system volume, S is the substrate sur- 
face area and p(r) the fluid density at point r. Con- 
sidering a simple van der Waals model [25], the bulk 
free energy per unit volume, V^b(p), can be given as 
Mp) = Pc{^P + - 2i/p + 3 - 2I3vt) in which 

= {P - Pc)/Pc and ut = {Tc - T)/Tc are the re- 
duced density and reduced temperature, respectively. 
The critical density, pressure, and temperature are set 
to pc — 7/2, pc = 1/8, and Tc = 4/7, respectively. Be- 
low Tc, the model describes liquid- vapor coexistence with 
related equilibrium densities pL,v = Pc(l i VP^t)- 

The parameter /3 is related to the interface thick ness 
^ and the surface tension a via ^ = ^J npl/ {A^vrpc) and 
CF = 4:/3^/2K,pc{f3lyT)^^'^pc' When combined with an ap- 
propriate variation of hz^ it allows to vary the surface 
tension and interface width independently. Using Cahn- 
Hilliard approach, the surface free energy per unit area, 
V^s, is approximated as — ^iPs, where ps is the density of 
fluid on the solid substrate and is a constant which 
can be used to tune the contact angle. Minimizing the 
free energy functional ^, Eq. ([2]), subject to the condition 
ips = — 0iPs leads to an equilibrium boundary condition 
for the spatial derivative of fluid density in the direction 
normal to the substrate, d±p = —(pi/n. The parameter 
is related to the Young contact angle via 



01 = 2/3rT\/2pc^sign(- - 6>)y cos - (1 - cos -), (3) 

where a = cos~^(sin^ ^y), and Oy is the equilibrium 
Young contact angle and the function "sign" determines 
the sign of its argument. All the quantities in this paper 



are given in dimensionless lattice Boltzmann units. The 
LB relaxation time is set to r = 0.8 and the tempera- 
ture is fixed to T = 0.4. For a typical choice of (5 = 0.1, 
for example, this leads to the equilibrium liquid and va- 
por densities pL ^ 4.1 and pv ^ 2.9. Depending on the 
case of interest, n lies in the range [0.002,0.008] and the 
size of the simulation box is varied with values around 
Lx X Ly X Lz = 125 X 90 X 90 lattice nodes for spher- 
ical, and Lx X Ly X Lz = 125 x 25 x 90 for cylindrical 
droplets. Periodic boundary conditions are applied in the 
X and ^-directions. 



III. DISCUSSION AND RESULTS 

As mentioned above, in the experimental reports which 
have considered the motion of a suspended droplet on 
surfaces with a gradient of texture, the behavior of 
droplets is not unique [H, [H, [2O, |2l| . This indicates that 
the roughness factor as well as the roughness density are 
not sufficient for a full characterization of a rough sur- 
face. Although Eq. ([1]) predicts a decrease of the effective 
contact angle upon an increase of roughness density and 
hence a driving force along the gradient of 0, the contact 
angle hysteresis [31] may be strong enough in order to 
prevent a spontaneous droplet motion [2l|. 

The present work underlines this aspect by explicitly 
showing that the behavior of a droplet on substrates pat- 
terned by pillar microstructure with the same pillar den- 
sity gradient, but different pillar geometries (e.g. rect- 
angular posts with different pillar width and spacing) 
can be qualitatively different. While in the one case the 
droplet spontaneously moves due to the roughness gradi- 
ent induced driving force, it may become arrested if an 
unfavorable geometry is chosen. 

A simple way to study the effect of a gradient tex- 
ture is to introduce an abrupt (stepwise) change in the 
roughness (pillar) density along a given spatial direction. 
Adopting this choice, we design a substrate divided into 
two regions, each with a constant pillar density. In order 
to underline the crucial role of pillar arrangement on the 
behavior of droplet, we consider two different cases of the 
same pillar density gradient as shown in Fig. [H 

Using the two substrates shown in Fig.[TJ we performed 
a series of lattice Boltzmann simulations placing at time 
t = a spherical liquid droplet close to the top of the 
border line separating the two regions of different pillar 
density. A close look at the left panels in Fig. [2] reveals 
that, both in the case of substrates A and B, the presence 
of a roughness gradient leads to an asymmetric spread- 
ing of droplet. However, despite this similarity of the 
dynamics at the early stages of spreading, the long time 
behavior of the droplet strongly depends on the specific 
arrangement of pillars. In particular, in the case of sub- 
strate A, the droplet motion is stopped on the gradient 
zone, while in case of substrate B it completely reaches 
the more favorable region of higher <p. 

In order to study the effect of droplet shape on the 



3 




FIG. 1: Top view of two step gradient substrates. In the left 
panel (referred to as case A) , the pillar density to the left side 
(x < 50) is (/)Left = 0.187 (square posts of length a = 6 = 3) 
while it is set to (/)Right = 0.321 on the right side (x > 50, 
rectangular posts of length a = 9 and width 6 = 3). The 
spacing distance of the pillars in the x-direction is = 5 
and in the ^-direction is = 3 overall on the substrate. The 
height of the posts is c = 6. The right panel (case B) is 
obtained from A by shifting the posts on each second raw 
horizontally by an amount of {a -\- dx)/2 with c^^, = 5, a = 3 
for X < 50 and a = 9 for x > 50. All lengths are given in LB 
units. 




FIG. 2: Initial setup and final states of a spherical droplet on 
substrates with an abrupt (step- wise) change of pillar density. 
The cases of substrates A and B (see Fig. [T]) are compared. 



above phenomenon, we also performed a set of simula- 
tions using a cylindrical droplet instead of a sphere. The 
results of these simulations, shown in Fig. O are in line 
with the case of spherical droplets. 

For further considerations, we use the substrate type 
B. The left panel of Fig. [H depicts the xz-cross section 
through the center of mass of a spherical drop {Rq^ = 40) 
at different times during its motion over the step gra- 
dient zone (the pillar densities on the left and right 
halves of the substrate are fixed to ^Left = 0.187 and 
^Right = 0.375). The interested reader can see the mo- 
tion of droplet in the supplementary movie. The corre- 
sponding footprint of the droplet (three phase contact 
line) is shown in the right panel of Fig. [H 

The footprint of the droplet reflects the geometry 




FIG. 3: The same set of simulations as in Fig. [2] but for the 
case of cylindrical droplets. 




20 40 60 80 100 120 20 40 60 80 100 



FIG. 4: The x^-cross section of the liquid-vapor interface 
(left) and the corresponding footprint (right) of a spherical 
droplet on a step gradient substrate. In the both panels, the 
time increases from left to right: t = 5 x lO'^, 6 x 10^, 1.2 x 10^ 
and 2 X 10^ 



(shape and arrangement) of the posts. A trend towards 
increasing contact area is also observed in accordance 
with the lower effective contact angle in the right region 
(higher pillar density). A closer look at the footprints 
in Fig. m (right panel) shows how the chess board-like 
arrangement of the posts allows the droplet to find the 
neighbor posts in the gradient direction. 

Next we create substrate patterns of type B with vari- 
ous values of by keeping ^Left unchanged and varying 
^Right- Results on the dynamics of a cylindrical drop on 
such texture gradient substrates are shown in Fig. [5l A 
survey of the center of mass position versus time in Fig. [5] 
reveals that the droplet motion is first linear in time un- 
til it reaches a constant value. The plateau corresponds 
to the case, where the droplet has completely left the re- 
gion of lower pillar density. Since no driving force exists 
in this state, the droplet velocity vanishes due to viscous 
dissipation. 

Using the linear part of the data shown in Fig. [5l we 
define an average velocity for the motion of droplet 's cen- 
ter of mass upon the action of texture gradient forces. 
Importantly, Fig. [5] reveals the strong effect of the sur- 
face tension on droplet dynamics. Both absolute values 
of droplet velocity for a given A(j) as well as the slope of 



4 



tn 
v> 
m 

E 
**— 
o 

0) 


o 



.80 
70 
60 
50 



# ,f|** R=40 

P - ' 

cylindrical droplet 



• s « 




\i 




■ A( 


t>=0.134 


> A( 


)=0.146 


• A( 


)=0.188 


A( 


)=0.163 



r]u/R'^ - (cr/i?2)V^A^. Hence, 



10 
time 



15 



X 10 



20 

5 



FIG. 5: The x-component of the center of mass position versus 
time for a cyUndrical droplet using, (/)Left = 0.187 and ^mght = 
0.321, 0.333, 0.35 and 0.375. The two group of curves belong 
to two different surface tensions of ctq = 5.4 x 10~^(LB units) 
(right; also labeled as (1) for further reference) and 4ao (left). 
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FIG. 6: Left: Droplet 's center of mass velocity versus the 
difference in pillar density, — (/)Right — ^Left , extracted from 
the linear part of the center of mass motion (see e.g. Fig. [5)). 
Results for three different liquid-vapor surface tensions are 
depicted. From top to bottom: 4cro, 2cro and ctq, where ctq = 
5.4 X 10~^ (LB units). In all cases, a linear variation is seen in 
accordance with the simple model, Eq. (|4|. Right: A further 
test of Eq. (|4|, where the dependence of droplet velocity on 
surface tension is shown for a fixed A(/). 



the data significantly depend upon a. 

In order to rationalize these observations, we provide 
a simple model based on scaling arguments. Noting that 
the flow we consider is in the viscous regime, we neglect 
inert ial terms in the Navier- Stokes equation and write 
for the steady state = — Vp + r/Aii, where u is the fluid 
velocity, p is the hydrostatic pressure and r] the viscosity. 
The velocity u varies only over a distance of the order 
of the droplet radius, hence A^x ~ ujB?. On the other 
hand, Vp -c^PLapiace('^£)/^ = {cr/R^)dR/R, assuming 
that the driving force originates from the Laplace pres- 
sure variation (over a length of the order of R) within 
the droplet. For the case of a cylindrical droplet of 
unit length, the condition of constant droplet volume, 
Vt = R^[e*^ - sin(2l9£)/2], Eq. ^ and some algebra 
lead to dR/R ~ {tt — OQ)d(j) ~ \^d(j) (the relation 
TT — Og^ ^ follows from Eq. ([T]) assuming Oq close 
to TT [2l|). Putting all together, and after a change 
of notation d(p = Acj) = (/)Right — (/^Left, we arrive at 



u ~ — v/^A0. 
V 



(4) 



Interestingly, despite different mechanisms at work, 
both Eq. (|4]) and Eq. (5) in [21] predict a linear depen- 
dence of droplet velocity on A^. In [2l| , the lateral veloc- 
ity is estimated from roughness gradient induced asym- 
metry of dewetting of a droplet, flattened due to impact. 
The situation we consider is different. There is no impact 
and hence a related flattening is absent in the present 
case. Furthermore, the dynamics we study is in the vis- 
cous regime whereas the high impact velocity in [21] sup- 
ports the relevance of inertia. These differences show up 
in different predictions regarding the dependence of the 
droplet velocity on surface tension, fluid viscosity and 
density. While Eq. (5) in [21] predicts a dependence on 
the square root of a, Eq. (|4|) suggests that, in our case, 
a linear dependence on a should be expected. 

We therefore examine Eq. (|4]) not only with regard to 
the relation between droplet velocity u and difference in 
roughness density A^ (left panel in Fig. [6]) but also check 
how u changes upon a variation of the surface tension a 
for a fixed Acj). Results of this latter test are depicted in 
the right panel of Fig. [6l confirming the expected linear 
dependence of on a. It is noteworthy that a in the right 
panel of Fig. [6] varies roughly by a factor of 10 so that a 
square root dependence can definitely be ruled out. 

It is worth emphasizing that the above discussed lin- 
ear relation between droplet velocity u and difference in 
pillar density A(j) is expected to hold as long as pinning 
forces are weaker than texture gradient induced driving 
force. In this case, the specific details of pillar arrange- 
ment seem to modify the prefactors entering the scal- 
ing relation, Eq. (|4]), but not the predicted linear law. 
Figure [3 is devoted to this aspect. In this figure, the 
dependence of u on A0 is compared for two slightly dif- 
ferent ways of realizing A^: In (1) 0Left = 0.187 while 
^Left = 0.2 in (2). All other aspects/parameters are iden- 
tical. As a consequence, since the list of investigated 
</>Right is exactly the same, slightly higher values of A<p 
are realized in (1) as compared to (2). If the details of pil- 
lar arrangement were unimportant, all the velocity data 
obtained from these two series of simulations should lie 
on the same line. As shown in Fig. [71 this is obviously 
not the case. Rather, the linear relation between u and 
A^ seems to hold independently for each studied case. 

Next we examine how the droplet 's contact area 
changes with time as the droplet moves on the gradient 
zone. We determine this quantity by simply counting 
the number of grid points beneath the droplet. The time 
evolution of the contact area is compared to that of the 
center of mass position in Fig. [8] (left panel) for a spheri- 
cal droplet of radius Rq = 36. In contrast to the center of 
mass position, which increases monotonously with time, 
the area beneath the droplet exhibits irregularities and 
oscillations. 

We presume that these irregularities and oscillations 
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FIG. 7: Left: The x-component of the center of mass position 
versus time for a cyhndrical droplet using, (/)Left = 0.2 and 
0Right = 0.321, 0.333, 0.35 and 0.375 [a = ctq = 5.4x10"^ (LB 
units)]. Right: Droplet 's center of mass velocity extracted 
from the linear part of the data shown in the left panel and 
in Fig. E] (labeled as (1)). 
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FIG. 8: Left: The contact area and the x-component of 
droplet 's center of mass position versus time. Right: Foot- 
print of the droplet on the step gradient substrate for times 
^A, and tc corresponding to the extrema of the contact 
area as labeled in the left panel. 



are closely related to the dynamics of three phase con- 
tact line. This idea is confirmed by the plot in the right 
panel of Fig. [51 where droplet 's footprint is shown for 
times ^A, and tc corresponding to the three extrema 
in the contact area labeled by A, B and C. As seen from 
this plot, the increase of the droplet 's base area between 
times ^A and is accompanied by a significant motion 
of the three phase contact line on the right side of the 
footprint while it remains essentially pinned to the pil- 
lars on the left side (with the exception of depinning from 
the left most pillar). The state B is, however, energeti- 
cally unfavorable due to a stretched shape of the contact 
line. The transition from B to C reduces this asymme- 
try, thereby leading to a smaller contact area at tc- We 
emphasize here that such local events are not included in 
the Cassie-Baxter picture, Eq. ([1]). 



IV. CONCLUSION 

We use a two-phase lattice Boltzmann model to study 
the dynamic behavior of suspended droplets on patterned 



hydrophobic substrates with a step-wise change in pil- 
lar density. We show that the specific arrangement of 
pillars may play a significant role for the dynamics of 
the droplet on such substrates. In particular, varying 
the pillar arrangement while keeping the gradient of pil- 
lar density unchanged (Fig. [1]), we show that both full 
transport over the gradient zone as well as complete ar- 
rest between the two regions of different pillar density 
may occur (Figs. [2] and [3]). 

The relation between the droplet motion and the gra- 
dient of pillar density is investigated, revealing a lin- 
ear dependence for the range of parameters studied 
(Figs. [5] and [6]). A simple model is provided based on the 
balance between the viscous dissipation and the driving 
force, the latter assumed as the gradient of the inter- 
nal droplet (Laplace) pressure (Eq. (j3|)). The model not 
only reproduces the observed dependence on the pillar 
density gradient but also predicts a linear dependence of 
the steady state droplet velocity on the surface tension. 
This prediction is in line with results of lattice Boltz- 
mann simulations, where the surface tension is varied by 
roughly a factor of 10 (Fig. [6|). 

Moreover, comparing droplet dynamics for two slightly 
different ways of realizing the gradient of texture, it 
is shown that the gradient in pillar density does not 
uniquely determine the droplet velocity. Rather, the way 
this gradient is implemented also matters to some extent 
(Fig. in. 

A detailed survey of the contact line dynamics is 
also provided revealing interesting pinning and depin- 
ning events leading to small amplitude oscillations of the 
droplet 's contact area during its motion over the gradient 
zone (Fig. [8|). 
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